clear all
set more off

// Robustness population (mergers, epci, large mkt)

global input "...\replication\data"
global output "...\replication\outputs"

use "$input\com_data.dta", clear

ren com_res com
collapse (firstnm) dep_res (mean) dist_border (sum) ind, by(com)

****************** Costmetics ******************
** Remove Corsica
drop if substr(com,2,1) == "A"
drop if substr(com,2,1) == "B"
destring com, replace


****************** Get dist border ******************
merge 1:1 com using "$input\com_arcgis.dta"
drop if _merge != 3
drop _merge
keep if dist_nat_border < 0 
keep if dist_river < 0


merge m:1 com using "$input\com_mergers.dta"
drop if _merge == 2
drop _merge
replace mergers = 0 if mergers == .


merge m:1 com using "$input\com_epci.dta"
drop if _merge == 2
drop _merge
replace epci = 0 if epci == .


****************** Get neighbor region for all com ******************
tostring com, replace
replace com = "0"+com if strlen(com) == 4
merge 1:1 com using "$input\loc_i_raw.dta"
drop if _merge == 2
drop _merge
destring com dep_res, replace

********* Other side ****************
geonear com x_com y_com using "$input\loc_j_raw.dta" , neighbors(dep x_dep y_dep) nearcount(2)
gen other_dep = nid1 if nid1 != dep_res
replace other_dep = nid2 if nid1 == dep_res

egen home_lab_mkt = sum(ind), by(dep_res)
egen other_lab_mkt = sum(ind), by(other_dep)


****************** Gen treatment ******************

replace dist_dep = -dist_dep if other_lab_mkt > home_lab_mkt

gen treat = (other_lab_mkt < home_lab_mkt)


//////////////////////////////////////////////////////////
//////////////// Excluding mergers ///////////////////////

****************** calonico et al ******************
rdbwselect ind dist_dep if mergers == 0, c(0) 
local opt_h=e(h_mserd)


rdrobust ind  dist_dep if mergers == 0,c(0) h(`opt_h') 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < `opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_1

rdrobust ind  dist_dep if mergers == 0,c(0) h(`opt_h')  covs(dep_res other_dep)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < `opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_2

rdrobust ind  dist_dep if mergers == 0,c(0) h(2*`opt_h')  covs(dep_res other_dep)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < 2*`opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_3

rdrobust ind  dist_dep if mergers == 0,c(0) h(0.5*`opt_h')  covs(dep_res other_dep)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < 0.5*`opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_4


****************** Results in Table ******************
esttab ATE_1 ATE_2 ATE_3 ATE_4, tex ///
b(%9.3f) se(%9.3f) stats( h_l N N_bw m_outc, labels("Bandwidth" "Obs" "Obs. bw." "Mean dep. var. below" ) ///
fmt( %9.3gc %9.0gc %9.0gc %9.3gc)) nolines nogaps star(* 0.10 ** 0.05 *** 0.01) ///
keep(RD_Estimate) coef(RD\_Estimate  "ATE") 


//////////////////////////////////////////////////////////
//////////////// Excluding epci ///////////////////////

****************** calonico et al ******************
rdbwselect ind dist_dep if epci == 0, c(0) 
local opt_h=e(h_mserd)


rdrobust ind  dist_dep if epci == 0,c(0) h(`opt_h') 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < `opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_1

rdrobust ind  dist_dep if epci == 0,c(0) h(`opt_h')  covs(dep_res other_dep)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < `opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_2

rdrobust ind  dist_dep if epci == 0,c(0) h(2*`opt_h')  covs(dep_res other_dep)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < 2*`opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_3

rdrobust ind  dist_dep if epci == 0,c(0) h(0.5*`opt_h')  covs(dep_res other_dep)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < 0.5*`opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_4


****************** Results in Table ******************
esttab ATE_1 ATE_2 ATE_3 ATE_4, tex ///
b(%9.3f) se(%9.3f) stats( h_l N N_bw m_outc, labels("Bandwidth" "Obs" "Obs. bw." "Mean dep. var. below" ) ///
fmt( %9.3gc %9.0gc %9.0gc %9.3gc)) nolines nogaps star(* 0.10 ** 0.05 *** 0.01) ///
keep(RD_Estimate) coef(RD\_Estimate  "ATE") 


//////////////////////////////////////////////////////////
//////////////// Excluding large labor markets ///////////////////////

gen large_mkt = 0
replace large_mkt = 1 if dep_res == 75
replace large_mkt = 1 if other_dep == 75
replace large_mkt = 1 if dep_res == 69
replace large_mkt = 1 if other_dep == 69
replace large_mkt = 1 if dep_res == 13
replace large_mkt = 1 if other_dep == 13

****************** calonico et al ******************
rdbwselect ind dist_dep if large_mkt == 0, c(0) 
local opt_h=e(h_mserd)


rdrobust ind  dist_dep if large_mkt == 0,c(0) h(`opt_h') 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < `opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_1

rdrobust ind  dist_dep if large_mkt == 0,c(0) h(`opt_h')  covs(dep_res other_dep)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < `opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_2

rdrobust ind  dist_dep if large_mkt == 0,c(0) h(2*`opt_h')  covs(dep_res other_dep)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < 2*`opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_3

rdrobust ind  dist_dep if large_mkt == 0,c(0) h(0.5*`opt_h')  covs(dep_res other_dep)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < 0.5*`opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_4


****************** Results in Table ******************
esttab ATE_1 ATE_2 ATE_3 ATE_4, tex ///
b(%9.3f) se(%9.3f) stats( h_l N N_bw m_outc, labels("Bandwidth" "Obs" "Obs. bw." "Mean dep. var. below" ) ///
fmt( %9.3gc %9.0gc %9.0gc %9.3gc)) nolines nogaps star(* 0.10 ** 0.05 *** 0.01) ///
keep(RD_Estimate) coef(RD\_Estimate  "ATE") 